#include "EdsCalcTask.h"
#include "mfAddFuncCalcTask.h"
#include "kernel/feObserver.h"
#include <math.h>
#include "Timer.h"
#include "FastMath.h"

//---------------------------------------------------------------------------------------------------------------------------------------------
inline double gCalcPhi(double iCenterX, double iCenterY, double iX, double iY)
{
	//get phi-coord in polar coords
	iX -= iCenterX;
	iY -= iCenterY;

	if(iX == 0 && iY > 0) return PI/2;
	if(iX == 0 && iY < 0) return 3*PI/2;

	double phi = atan(fabs(iY/iX));
	if(iX >= 0 && iY >= 0) return phi;
	if(iX <= 0 && iY >= 0) return PI - phi;
	if(iX >= 0 && iY <= 0) return 2*PI - phi;	
	if(iX <= 0 && iY <= 0) return phi + PI;

	return phi;
}


//---------------------------------------------------------------------------------------------------------------------------------------------
inline double gCalcR(double iX1, double iY1, double iX2, double iY2)
{
	//calc distance btw 2 points
	iX2 -= iX1;
	iY2 -= iY1;
	return FastMath::Sqrt(iX2*iX2 + iY2*iY2);
}

//---------------------------------------------------------------------------------------------------------------------------------------------
VmdEdsCalcTask::VmdEdsCalcTask()
{
	Layers = 0;
	NumLayers = 0;
	TimesArray = 0;
	TimesCount = 0;
}

//---------------------------------------------------------------------------------------------------------------------------------------------
double J0 (double iX) 
{
	if (iX == 0) return 1.0;

	const double
		P1=1.0, P2=-0.1098628627E-2, P3=0.2734510407E-4,
		P4=-0.2073370639E-5, P5= 0.2093887211E-6,
		Q1=-0.1562499995E-1, Q2= 0.1430488765E-3, Q3=-0.6911147651E-5,
		Q4= 0.7621095161E-6, Q5=-0.9349451520E-7,
		R1= 57568490574.0, R2=-13362590354.0, R3=651619640.7,
		R4=-11214424.18, R5= 77392.33017, R6=-184.9052456,
		S1= 57568490411.0, S2=1029532985.0, S3=9494680.718,
		S4= 59272.64853, S5=267.8532712, S6=1.0;
	double
		AX,FR,FS,Z,FP,FQ,XX,Y, TMP;


	AX = fabs(iX);
	if (AX < 8.0) 
	{
		Y = iX*iX;
		FR = R1+Y*(R2+Y*(R3+Y*(R4+Y*(R5+Y*R6))));
		FS = S1+Y*(S2+Y*(S3+Y*(S4+Y*(S5+Y*S6))));
		TMP = FR/FS;
	}
	else 
	{
		Z = 8./AX;
		Y = Z*Z;
		XX = AX-0.785398164;
		FP = P1+Y*(P2+Y*(P3+Y*(P4+Y*P5)));
		FQ = Q1+Y*(Q2+Y*(Q3+Y*(Q4+Y*Q5)));

		double s;
		double c;

		FastMath::SinCos(XX, s, c);		
		TMP = FastMath::Sqrt(0.636619772/AX)*(FP*c-Z*FQ*s);
	}
	return TMP;
}

//---------------------------------------------------------------------------------------------------------------------------------------------
bool VmdEdsCalcTask::Init(ConductiveLayer *iLayers, int iNumLayers, fePolyline<double> *iGenerator, fePolyline<double> *iReceiver)
{
	if(Layers)
		delete [] Layers;

	//create and copy layers
	NumLayers = iNumLayers;
	Layers = new ConductiveLayer[iNumLayers];
	if(!Layers) return false;
	for(int i=0; i<NumLayers; i++)
		Layers[i] = iLayers[i];
	
	Generator = iGenerator;
	Receiver = iReceiver;


	//create triangulation
	RecTris.Init(Receiver);
	RecTris.Triangulate();
	RecTris.Tesselate(5);
	RecTris.FillAreas();

	Dipole.SetMedium(Layers, NumLayers);

	return true;
}


//---------------------------------------------------------------------------------------------------------------------------------------------
double VmdEdsCalcTask::Integrate(double iF[3], double iX[3])
{
/*
	double t1 = iX[2]*iX[2];
	double t2 = t1*iX[2];
	double t14 = iF[0]*t1;
	double t18 = iX[2]*iF[0];
	double t19 = iX[0]*iX[0];
	double t22 = iX[1]*iX[1];
	double t37 = iF[2]*iX[1];
	double t45 = t19*iX[0];
	double t54 = iF[0]*t2-iF[1]*t2+3.0*t1*iF[1]*iX[0]-2.0*t1*iF[2]*iX[1]+2.0*iX[0]*iF[2]*t1
		-4.0*t14*iX[1]+t14*iX[0]-2.0*t18*t19+3.0*t18*t22-iF[2]*t19*iX[2]+2.0*t18*iX[1]*iX[0]+
		3.0*iX[2]*iF[2]*t22-3.0*iF[1]*iX[2]*t19-2.0*iX[2]*iX[0]*t37+4.0*t37*t19-3.0*iF[2]*t22*
		iX[0]+iF[1]*t45-3.0*iF[0]*t22*iX[0]-iF[2]*t45+2.0*iF[0]*iX[1]*t19;
	return t54/(iX[1]-iX[0])/(iX[1]-iX[2])/6.0;
*/


	return 0.5*(iF[0]+iF[1])*(iX[1]-iX[0]);
}


//---------------------------------------------------------------------------------------------------------------------------------------------
bool VmdEdsCalcTask::Eds(double iZ, double *iTimes, int iTimesCount, double *oResV)
{
	//check layers
	if(NumLayers < 1 || !Layers || !iTimes || iTimesCount < 1) return false;

	//fill time data
	TimesArray = iTimes;
	TimesCount = iTimesCount;

	ZeroMemory(oResV, sizeof(oResV[0])*iTimesCount);


	double *dbzdt1 = new double[TimesCount];
	double *dbzdt2 = new double[TimesCount];
	double *dbzdt3 = new double[TimesCount];


	for(int trinum=0; trinum<RecTris.GetTriNum(); trinum++)
	{
		//get tri
		fePolyline<double>::V2 tri[3];
		RecTris.GetTri(trinum, tri[0], tri[1], tri[2]);


		//integrate eds
		fePolyline<double>::V2 p;
		p.X = (tri[0].X + tri[1].X + tri[2].X) / 3;
		p.Y = (tri[0].Y + tri[1].Y + tri[2].Y) / 3;


		//get bbox in (phi, r) coords
		double minr = 1e40;
		double maxr = 0;
		double minphi = 20;
		double maxphi = 0;
		for(int i=0; i<Generator->GetSize(); i++)
		{
			fePolyline<double>::V2 p1 = Generator->GetPoint(i);
			double phi = gCalcPhi(p.X, p.Y, p1.X, p1.Y);
			double r = gCalcR(p.X, p.Y, p1.X, p1.Y);

			if(minphi > phi) minphi = phi;
			if(maxphi < phi) maxphi = phi;
			if(minr > r) minr = r;
			if(maxr < r) maxr = r;
		}
		
		if(Generator->TestPoint(&p)) 
		{
			minr = 0;
			minphi = 0;
			maxphi = 2*PI;
		}

		//get steps
		double dr = min(maxr-minr, 0.1);
		double dphi = min(maxphi-minphi, 0.1)/(maxr);
		int numrings = int((maxr-minr)/dr+0.5);
		int numphi = int((maxphi-minphi)/dphi+0.5);
		if(numphi < 500) 
		{
			dphi = (maxphi-minphi)/500;
			numphi = 500;
		}


		//get first dbzdt
		Dipole.DBzDt(minr == 0 ? 1e-5 : minr, iZ, dbzdt1, TimesArray, TimesCount);


		//integrate 
		for(int i=0; i<numrings; i++)
		{
			//get generator area in r..r+dr
			double r = minr + dr * i;
			double s = (maxphi-minphi) * ((r+dr)*(r+dr) - r*r)/(2*numphi);
			int indicator = 0;
			for(int nphi=0; nphi < numphi; nphi++)
			{
				double phi = minphi + dphi * nphi + 0.5*dphi;
				fePolyline<double>::V2 point;
				
				double s;
				double c;
				FastMath::SinCos(phi, s, c);
				point.X = (r + dr*0.5)*c + p.X;
				point.Y = (r + dr*0.5)*s + p.Y;
				if(Generator->TestPoint(&point)) indicator++;
			}
			s *= indicator/dr; //final area of generator


			//for simpson need 3 points
			Dipole.DBzDt(r+dr/2, iZ, dbzdt2, TimesArray, TimesCount);
			Dipole.DBzDt(r+dr, iZ, dbzdt3, TimesArray, TimesCount);

			double c = s * RecTris.GetArea(trinum);
			for(int it=0; it<TimesCount; it++)
			{
				//simpson
				oResV[it] += dr/6 * c * (dbzdt1[it] + 4 * dbzdt2[it] + dbzdt3[it]);
			}

			double *b = dbzdt1;
			dbzdt1 = dbzdt3;
			dbzdt3 = b;
		}
	}


//	fePolylineGenerator g;
//	g.SavePolylineToFile((fePolyline<float> *)&eds, "eds_simple");
	
//	feObserverManager::GetInstance().Notify("Time to solve: %f sec", t.Seconds());

	delete [] dbzdt1;
	delete [] dbzdt2;
	delete [] dbzdt3;

	return true;
}


//---------------------------------------------------------------------------------------------------------------------------------------------
VmdEdsCalcTask::~VmdEdsCalcTask()
{

}
